home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Collection of Tools & Utilities
/
Collection of Tools and Utilities.iso
/
fortran
/
linpklib.zip
/
SSVDC.FOR
< prev
next >
Wrap
Text File
|
1984-01-06
|
15KB
|
481 lines
SUBROUTINE SSVDC(X,LDX,N,P,S,E,U,LDU,V,LDV,WORK,JOB,INFO)
INTEGER LDX,N,P,LDU,LDV,JOB,INFO
REAL X(LDX,1),S(1),E(1),U(LDU,1),V(LDV,1),WORK(1)
C
C
C SSVDC IS A SUBROUTINE TO REDUCE A REAL NXP MATRIX X BY
C ORTHOGONAL TRANSFORMATIONS U AND V TO DIAGONAL FORM. THE
C DIAGONAL ELEMENTS S(I) ARE THE SINGULAR VALUES OF X. THE
C COLUMNS OF U ARE THE CORRESPONDING LEFT SINGULAR VECTORS,
C AND THE COLUMNS OF V THE RIGHT SINGULAR VECTORS.
C
C ON ENTRY
C
C X REAL(LDX,P), WHERE LDX.GE.N.
C X CONTAINS THE MATRIX WHOSE SINGULAR VALUE
C DECOMPOSITION IS TO BE COMPUTED. X IS
C DESTROYED BY SSVDC.
C
C LDX INTEGER.
C LDX IS THE LEADING DIMENSION OF THE ARRAY X.
C
C N INTEGER.
C N IS THE NUMBER OF COLUMNS OF THE MATRIX X.
C
C P INTEGER.
C P IS THE NUMBER OF ROWS OF THE MATRIX X.
C
C LDU INTEGER.
C LDU IS THE LEADING DIMENSION OF THE ARRAY U.
C (SEE BELOW).
C
C LDV INTEGER.
C LDV IS THE LEADING DIMENSION OF THE ARRAY V.
C (SEE BELOW).
C
C WORK REAL(N).
C WORK IS A SCRATCH ARRAY.
C
C JOB INTEGER.
C JOB CONTROLS THE COMPUTATION OF THE SINGULAR
C VECTORS. IT HAS THE DECIMAL EXPANSION AB
C WITH THE FOLLOWING MEANING
C
C A.EQ.0 DO NOT COMPUTE THE LEFT SINGULAR
C VECTORS.
C A.EQ.1 RETURN THE N LEFT SINGULAR VECTORS
C IN U.
C A.GE.2 RETURN THE FIRST MIN(N,P) SINGULAR
C VECTORS IN U.
C B.EQ.0 DO NOT COMPUTE THE RIGHT SINGULAR
C VECTORS.
C B.EQ.1 RETURN THE RIGHT SINGULAR VECTORS
C IN V.
C
C ON RETURN
C
C S REAL(MM), WHERE MM=MIN(N+1,P).
C THE FIRST MIN(N,P) ENTRIES OF S CONTAIN THE
C SINGULAR VALUES OF X ARRANGED IN DESCENDING
C ORDER OF MAGNITUDE.
C
C E REAL(P).
C E ORDINARILY CONTAINS ZEROS. HOWEVER SEE THE
C DISCUSSION OF INFO FOR EXCEPTIONS.
C
C U REAL(LDU,K), WHERE LDU.GE.N. IF JOBA.EQ.1 THEN
C K.EQ.N, IF JOBA.GE.2 THEN
C K.EQ.MIN(N,P).
C U CONTAINS THE MATRIX OF RIGHT SINGULAR VECTORS.
C U IS NOT REFERENCED IF JOBA.EQ.0. IF N.LE.P
C OR IF JOBA.EQ.2, THEN U MAY BE IDENTIFIED WITH X
C IN THE SUBROUTINE CALL.
C
C V REAL(LDV,P), WHERE LDV.GE.P.
C V CONTAINS THE MATRIX OF RIGHT SINGULAR VECTORS.
C V IS NOT REFERENCED IF JOB.EQ.0. IF P.LE.N,
C THEN V MAY BE IDENTIFIED WITH X IN THE
C SUBROUTINE CALL.
C
C INFO INTEGER.
C THE SINGULAR VALUES (AND THEIR CORRESPONDING
C SINGULAR VECTORS) S(INFO+1),S(INFO+2),...,S(M)
C ARE CORRECT (HERE M=MIN(N,P)). THUS IF
C INFO.EQ.0, ALL THE SINGULAR VALUES AND THEIR
C VECTORS ARE CORRECT. IN ANY EVENT, THE MATRIX
C B = TRANS(U)*X*V IS THE BIDIAGONAL MATRIX
C WITH THE ELEMENTS OF S ON ITS DIAGONAL AND THE
C ELEMENTS OF E ON ITS SUPER-DIAGONAL (TRANS(U)
C IS THE TRANSPOSE OF U). THUS THE SINGULAR
C VALUES OF X AND B ARE THE SAME.
C
C LINPACK. THIS VERSION DATED 03/19/79 .
C G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.
C
C ***** USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS.
C
C EXTERNAL SROT
C BLAS SAXPY,SDOT,SSCAL,SSWAP,SNRM2,SROTG
C FORTRAN ABS,AMAX1,MAX0,MIN0,MOD,SQRT
C
C INTERNAL VARIABLES
C
INTEGER I,ITER,J,JOBU,K,KASE,KK,L,LL,LLS,LM1,LP1,LS,LU,M,MAXIT,
* MM,MM1,MP1,NCT,NCTP1,NCU,NRT,NRTP1
REAL SDOT,T,R
REAL B,C,CS,EL,EMM1,F,G,SNRM2,SCALE,SHIFT,SL,SM,SN,SMM1,T1,TEST,
* ZTEST
LOGICAL WANTU,WANTV
C
C
C SET THE MAXIMUM NUMBER OF ITERATIONS.
C
MAXIT = 30
C
C DETERMINE WHAT IS TO BE COMPUTED.
C
WANTU = .FALSE.
WANTV = .FALSE.
JOBU = MOD(JOB,100)/10
NCU = N
IF (JOBU .GT. 1) NCU = MIN0(N,P)
IF (JOBU .NE. 0) WANTU = .TRUE.
IF (MOD(JOB,10) .NE. 0) WANTV = .TRUE.
C
C REDUCE X TO BIDIAGONAL FORM, STORING THE DIAGONAL ELEMENTS
C IN S AND THE SUPER-DIAGONAL ELEMENTS IN E.
C
INFO = 0
NCT = MIN0(N-1,P)
NRT = MAX0(0,MIN0(P-2,N))
LU = MAX0(NCT,NRT)
IF (LU .LT. 1) GO TO 170
DO 160 L = 1, LU
LP1 = L + 1
IF (L .GT. NCT) GO TO 20
C
C COMPUTE THE TRANSFORMATION FOR THE L-TH COLUMN AND
C PLACE THE L-TH DIAGONAL IN S(L).
C
S(L) = SNRM2(N-L+1,X(L,L),1)
IF (S(L) .EQ. 0.0E0) GO TO 10
IF (X(L,L) .NE. 0.0E0) S(L) = SIGN(S(L),X(L,L))
CALL SSCAL(N-L+1,1.0E0/S(L),X(L,L),1)
X(L,L) = 1.0E0 + X(L,L)
10 CONTINUE
S(L) = -S(L)
20 CONTINUE
IF (P .LT. LP1) GO TO 50
DO 40 J = LP1, P
IF (L .GT. NCT) GO TO 30
IF (S(L) .EQ. 0.0E0) GO TO 30
C
C APPLY THE TRANSFORMATION.
C
T = -SDOT(N-L+1,X(L,L),1,X(L,J),1)/X(L,L)
CALL SAXPY(N-L+1,T,X(L,L),1,X(L,J),1)
30 CONTINUE
C
C PLACE THE L-TH ROW OF X INTO E FOR THE
C SUBSEQUENT CALCULATION OF THE ROW TRANSFORMATION.
C
E(J) = X(L,J)
40 CONTINUE
50 CONTINUE
IF (.NOT.WANTU .OR. L .GT. NCT) GO TO 70
C
C PLACE THE TRANSFORMATION IN U FOR SUBSEQUENT BACK
C MULTIPLICATION.
C
DO 60 I = L, N
U(I,L) = X(I,L)
60 CONTINUE
70 CONTINUE
IF (L .GT. NRT) GO TO 150
C
C COMPUTE THE L-TH ROW TRANSFORMATION AND PLACE THE
C L-TH SUPER-DIAGONAL IN E(L).
C
E(L) = SNRM2(P-L,E(LP1),1)
IF (E(L) .EQ. 0.0E0) GO TO 80
IF (E(LP1) .NE. 0.0E0) E(L) = SIGN(E(L),E(LP1))
CALL SSCAL(P-L,1.0E0/E(L),E(LP1),1)
E(LP1) = 1.0E0 + E(LP1)
80 CONTINUE
E(L) = -E(L)
IF (LP1 .GT. N .OR. E(L) .EQ. 0.0E0) GO TO 120
C
C APPLY THE TRANSFORMATION.
C
DO 90 I = LP1, N
WORK(I) = 0.0E0
90 CONTINUE
DO 100 J = LP1, P
CALL SAXPY(N-L,E(J),X(LP1,J),1,WORK(LP1),1)
100 CONTINUE
DO 110 J = LP1, P
CALL SAXPY(N-L,-E(J)/E(LP1),WORK(LP1),1,X(LP1,J),1)
110 CONTINUE
120 CONTINUE
IF (.NOT.WANTV) GO TO 140
C
C PLACE THE TRANSFORMATION IN V FOR SUBSEQUENT
C BACK MULTIPLICATION.
C
DO 130 I = LP1, P
V(I,L) = E(I)
130 CONTINUE
140 CONTINUE
150 CONTINUE
160 CONTINUE
170 CONTINUE
C
C SET UP THE FINAL BIDIAGONAL MATRIX OR ORDER M.
C
M = MIN0(P,N+1)
NCTP1 = NCT + 1
NRTP1 = NRT + 1
IF (NCT .LT. P) S(NCTP1) = X(NCTP1,NCTP1)
IF (N .LT. M) S(M) = 0.0E0
IF (NRTP1 .LT. M) E(NRTP1) = X(NRTP1,M)
E(M) = 0.0E0
C
C IF REQUIRED, GENERATE U.
C
IF (.NOT.WANTU) GO TO 300
IF (NCU .LT. NCTP1) GO TO 200
DO 190 J = NCTP1, NCU
DO 180 I = 1, N
U(I,J) = 0.0E0
180 CONTINUE
U(J,J) = 1.0E0
190 CONTINUE
200 CONTINUE
IF (NCT